SoSe2021
\[Weight_{ij}=\alpha_i+\beta_i*Age_j + \epsilon_{ij}\]
\[Weight_{female,j}=\alpha_{female}+\beta_{female}*Age_j + \epsilon_{female,j}\] \[Weight_{male,j}=\alpha_{male}+\beta_{male}*Age_j + \epsilon_{male,j}\]
ChinookArgLängen und Gewichte für Königslachse von drei Orten in Argentinien (Daten sind im FSA Paket enthalten):
data("ChinookArg", package = "FSA")
str(ChinookArg)
'data.frame': 112 obs. of 3 variables: $ tl : num 120 115 111 110 110 ... $ w : num 17.9 17.2 16.8 15.8 14.3 13.8 12.8 11.7 12.8 14.8 ... $ loc: Factor w/ 3 levels "Argentina","Petrohue",..: 1 1 1 1 1 1 1 1 1 1 ...
Da für die meisten Arten die Beziehung zwischen dem Gewicht eines Tieres und seiner Länge einer 2-Parameter Powerfunktion (\(W = aL^{b}\)) folgt, wollen wir die Variablen log-transformieren:
ChinookArg$tl_log <- log(ChinookArg$tl) ChinookArg$w_log <- log(ChinookArg$w)
full_mod <- lm(w_log ~ loc + tl_log + loc:tl_log, ChinookArg) # oder: ~ loc * tl_log summary(full_mod)
Call:
lm(formula = w_log ~ loc + tl_log + loc:tl_log, data = ChinookArg)
Residuals:
Min 1Q Median 3Q Max
-0.58273 -0.18471 -0.00186 0.13088 1.63620
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.6750 1.5904 -4.197 5.64e-05 ***
locPetrohue -2.3957 3.1494 -0.761 0.449
locPuyehue -2.0696 1.6868 -1.227 0.223
tl_log 1.9836 0.3530 5.619 1.56e-07 ***
locPetrohue:tl_log 0.4795 0.6928 0.692 0.490
locPuyehue:tl_log 0.3624 0.3793 0.955 0.342
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3201 on 106 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8924
F-statistic: 185 on 5 and 106 DF, p-value: < 2.2e-16Falls x kategorial ist, macht es nicht viel Sinn das Model \(y = a + bx\) zu benutzen, da \(x\) nicht mit \(b\) multipliziert werden kann. Aber R hat eine Hilfskonstruktion:
model.matrix()Die Daten
df
y length 1 6.7 M 2 16.6 M 3 12.8 M 4 16.1 M 5 4.2 L 6 1.3 S 7 18.6 L 8 5.0 S 9 8.4 S 10 2.0 S
Die Dummy-Variablen
modelr::model_matrix(df, y ~ length)
# A tibble: 10 x 3
`(Intercept)` lengthM lengthS
<dbl> <dbl> <dbl>
1 1 1 0
2 1 1 0
3 1 1 0
4 1 1 0
5 1 0 0
6 1 0 1
7 1 0 0
8 1 0 1
9 1 0 1
10 1 0 1Wo ist die Längen-Klasse L geblieben?
model.matrix()Die Daten
df
y length 1 6.7 M 2 16.6 M 3 12.8 M 4 16.1 M 5 4.2 L 6 1.3 S 7 18.6 L 8 5.0 S 9 8.4 S 10 2.0 S
Die Dummy-Variablen
modelr::model_matrix(df, y ~ length)
# A tibble: 10 x 3
`(Intercept)` lengthM lengthS
<dbl> <dbl> <dbl>
1 1 1 0
2 1 1 0
3 1 1 0
4 1 1 0
5 1 0 0
6 1 0 1
7 1 0 0
8 1 0 1
9 1 0 1
10 1 0 1L ist durch den Achsenschnittpunkt (erste Spalte) dargestellt!
model.matrix(w_log ~ loc, ChinookArg) %>%
as.data.frame(.) %>%
# Hinzufügen der Original loc Variable
mutate(sampl_loc = ChinookArg$loc) %>%
head()
(Intercept) locPetrohue locPuyehue sampl_loc 1 1 0 0 Argentina 2 1 0 0 Argentina 3 1 0 0 Argentina 4 1 0 0 Argentina 5 1 0 0 Argentina 6 1 0 0 Argentina
Wir können sehen, dass die ersten 5 Lachse in ‘Argentinien’ gefangen wurden und dass die Beobachtungen mit einer 1 für die ‘intercept’ Variable (erste Faktorstufe → hier Argentinien) kodiert und 0 für die Variable locPetrohue und locPuyehue wurden.
Wie würdest Du diese Koeffizienten interpretieren (wenn also nur der Fangort die erklärende Variable ist)?
chin_mod <- lm(w_log ~ loc, ChinookArg) coef(chin_mod)
(Intercept) locPetrohue locPuyehue 2.25605528 -0.09844412 -1.52993775
Wie lautet die Gleichung für den Fangort?
Setze die berechneten Koeffizienten in die Regressionsgleichung mit Dummy-Variablen ein:
summary(full_mod)
Call:
lm(formula = w_log ~ loc + tl_log + loc:tl_log, data = ChinookArg)
Residuals:
Min 1Q Median 3Q Max
-0.58273 -0.18471 -0.00186 0.13088 1.63620
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.6750 1.5904 -4.197 5.64e-05 ***
locPetrohue -2.3957 3.1494 -0.761 0.449
locPuyehue -2.0696 1.6868 -1.227 0.223
tl_log 1.9836 0.3530 5.619 1.56e-07 ***
locPetrohue:tl_log 0.4795 0.6928 0.692 0.490
locPuyehue:tl_log 0.3624 0.3793 0.955 0.342
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3201 on 106 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8924
F-statistic: 185 on 5 and 106 DF, p-value: < 2.2e-16\[w.log_{ij} = \alpha_{i} + \beta_{i}*tl.log_{j} + \epsilon_{ij}\]
wobei i = Fangort, j = jeder einzelne Lachs
\[w.log_{A,j} = \alpha_{A} + \beta_{A}*tl.log_{j} + \epsilon_{A,j}\] \[w.log_{Pe,j} = \alpha_{Pe} + \beta_{Pe}*tl.log_{j} + \epsilon_{Pe,j}\] \[w.log_{Pu,j} = \alpha_{Pu} + \beta_{Pu}*tl.log_{j} + \epsilon_{Pu,j}\]
\[w.log = a_A + a_{Pe}*locPetrohue + a_{Pu}*locPuyehue +\\ b_A*tl.log+ b_{Pe}*locPetrohue*tl.log + b_{Pu}*locPuyehue*tl.log\]
coef(full_mod)
(Intercept) locPetrohue locPuyehue tl_log
-6.6749719 -2.3957331 -2.0695530 1.9836016
locPetrohue:tl_log locPuyehue:tl_log
0.4794543 0.3624014 Argentina: \(w.log_{A} = -6.67+1.98*tl_log\)
Petrohuye: \(w.log_{Pe}= (-6.67+-2.40)+(1.98+0.48*1)*tl.log = -9.07+2.46*tl.log\)
Puyehuye: \(w.log_{Pu} = (-6.67+-2.07)+(1.98+0.36*1)*tl.log = -8.74+2.34*tl.log\)
Anstelle der ‘all subsets selection’ (siehe Bsp. airquality), wollen wir hier etwas systematischer vorgehen und die sog. ‘backward selection’ vornehmen.
step() → führt eine vollständige Selektion durch, bis das Modell mit dem niedrigsten AIC gefunden wird.
method kann eingestellt werden, ob die Richtung ‘forward’, ‘backward’ oder beides (‘both’) sein soll.drop1 → erstellt reduzierte Modelle des übergebenen vollen Modells und gibt eine Vergleichstabelle mit den AICs zurück.par(mfrow = c(2,2)) plot(full_mod)
ChinookArg$residuals <- resid(full_mod) par(mfrow = c(1,2), mar = c(4,4,1,1)) boxplot(residuals ~loc, ChinookArg, main = "Residuals vs. loc"); abline(h = 0, lty = 2) plot(residuals ~ tl_log, ChinookArg, main = "Residuals vs. tl_log"); abline(h = 0, lty = 2)
Wenn Interaktionsterme im Modell vorkommen, müssen auch die enthaltenen Variablen als Hauptterme vorkommen. Es darf also nicht tl_log oder loc im reduzierten Modell fehlen. Nur die Interaktion selbst kann aktuell rausgenommen werden:
red_mod <- update(full_mod, .~.-loc:tl_log) AIC(full_mod, red_mod)
df AIC full_mod 7 70.53732 red_mod 5 67.57480
# Das reduzierte Modell wird das volle: full_mod <- red_mod # Einzelne X entfernen red_mod1 <- update(full_mod, .~.-loc) red_mod2 <- update(full_mod, .~.-tl_log) AIC(full_mod, red_mod1, red_mod2)
df AIC full_mod 5 67.57480 red_mod1 3 87.69126 red_mod2 4 224.06339
step(full_mod, method = "backward")
Start: AIC=-249.3
w_log ~ loc + tl_log + loc:tl_log
Df Sum of Sq RSS AIC
- loc:tl_log 2 0.1011 10.965 -252.27
<none> 10.864 -249.31
Step: AIC=-252.27
w_log ~ loc + tl_log
Df Sum of Sq RSS AIC
<none> 10.965 -252.267
- loc 2 2.634 13.599 -232.151
- tl_log 1 34.175 45.140 -95.779
Call:
lm(formula = w_log ~ loc + tl_log, data = ChinookArg)
Coefficients:
(Intercept) locPetrohue locPuyehue tl_log
-8.1217 -0.2281 -0.4570 2.3049 best_mod <- lm(w_log ~ loc + tl_log, ChinookArg) par(mfrow = c(2,2)) plot(best_mod)
ChinookArg$residuals <- resid(best_mod) par(mfrow = c(1,2), mar = c(4,4,1,1)) boxplot(residuals ~loc, ChinookArg, main = "Residuals vs. loc"); abline(h = 0, lty = 2) plot(residuals ~ tl_log, ChinookArg, main = "Residuals vs. tl_log"); abline(h = 0, lty = 2)
summary(best_mod)
Call:
lm(formula = w_log ~ loc + tl_log, data = ChinookArg)
Residuals:
Min 1Q Median 3Q Max
-0.60633 -0.19351 -0.00076 0.14352 1.61286
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.12169 0.56827 -14.292 < 2e-16 ***
locPetrohue -0.22813 0.08013 -2.847 0.00528 **
locPuyehue -0.45699 0.09231 -4.951 2.74e-06 ***
tl_log 2.30492 0.12563 18.347 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3186 on 108 degrees of freedom
Multiple R-squared: 0.8962, Adjusted R-squared: 0.8934
F-statistic: 311 on 3 and 108 DF, p-value: < 2.2e-16# Extraktion der Koeffizientenschätzung mit tidy() aus dem broom Paket broom::tidy(best_mod)
# A tibble: 4 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) -8.12 0.568 -14.3 1.19e-26 2 locPetrohue -0.228 0.0801 -2.85 5.28e- 3 3 locPuyehue -0.457 0.0923 -4.95 2.74e- 6 4 tl_log 2.30 0.126 18.3 5.72e-35
Um die ANOVA Tabelle zu erhalten, wo alle Faktorstufen miteinander verglichen werden, nutze die Funktion anova() für dein lm Objekt:
anova(best_mod)
Analysis of Variance Table
Response: w_log
Df Sum Sq Mean Sq F value Pr(>F)
loc 2 60.542 30.271 298.16 < 2.2e-16 ***
tl_log 1 34.175 34.175 336.61 < 2.2e-16 ***
Residuals 108 10.965 0.102
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1→ Dieser p-Wert für loc sollte im Bericht stehen, nicht die aus der summary Funktion!
coef(best_mod)
(Intercept) locPetrohue locPuyehue tl_log -8.1216871 -0.2281264 -0.4569901 2.3049203
Argentina: \(w.log_{A} = -8.12+2.30*tl.log\)
Petrohuye: \(w.log_{Pe}= (-8.12+-0.23)+2.30*tl.log = -8.35+2.30*tl.log\)
Puyehuye: \(w.log_{Pu} = (-8.12+-0.46)+2.30*tl.log = -8.58+2.30*tl.log\)
ChinookArg %>% data_grid(
tl_log = seq_range(tl_log, 20),
loc = factor(levels(loc))
) %>%
add_predictions(best_mod) %>%
ggplot(aes(x = tl_log, group = loc,
colour = loc)) +
geom_line(aes(y = pred)) +
geom_point(data = ChinookArg,
mapping = aes(y = w_log)) +
ylab("w_log") +
scale_colour_manual(values =
c("forestgreen", "orange1",
"deeppink4")) +
theme_classic()Bei weiteren Fragen: saskia.otto(at)uni-hamburg.de

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License except for the borrowed and mentioned with proper source: statements.
Image on title and end slide: Section of an infrared satellite image showing the Larsen C ice shelf on the Antarctic Peninsula - USGS/NASA Landsat: A Crack of Light in the Polar Dark, Landsat 8 - TIRS, June 17, 2017 (under CC0 license)